##################################################
# REPLICATION CODE
# "Unequal Policy Responsiveness in Spain"
# By Noam Lupu and Alejandro Tirado Castro
##################################################

rm (list = ls(all=TRUE))
library(foreign)
library(haven)
library(lme4)
library(ggplot2)
library(gridExtra)
library(ggeffects)
library(jtools)
library(robust)
library(performance)
library(interplot)
library(extrafont)


Dataset_CIS <- read_dta("CIS_data.dta")
Percentages <- read_dta("Demographics.dta")
attach (Dataset_CIS)


##################################################
#################### FIGURE A1 ###################
##################################################

jpeg(file="figurea1.jpg", width=114, height=80, units="mm", res=900, pointsize=7)

ggplot(Percentages, aes(x=year, y=primary)) + 
  geom_line(data = Percentages, aes(x=year, y=primary, col = "Primary", linetype = "Primary"), size=0.5) +
  geom_line(data = Percentages, aes(x=year, y=secondary, col = "Secondary", linetype="Secondary"), size=0.5) +
  geom_line(data = Percentages, aes(x=year, y=tertiary, col = "Tertiary", linetype="Tertiary"), size=0.5) +
  ylab("Proportion") + xlab("") + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.25)
  ) + ggtitle("") + theme(plot.title = element_text(hjust = 0.5)) + theme(text = element_text(size = 7, family="sans")) + 
  coord_cartesian(ylim = c(0,100)) +
  theme(legend.position="bottom") + theme(legend.key.width = unit(1,"cm")) +
  scale_color_manual(name = "", values = c("gray60", "gray30", "black"), labels = c("Primary", "Secondary", "Tertiary")) +
  scale_linetype_manual(name = "", values = c("solid", "dotted", "dotdash"), labels = c("Primary", "Secondary", "Tertiary")) 

dev.off()


##################################################
#################### FIGURE A2 ###################
##################################################

jpeg(file="figurea2.jpg", width=114, height=80, units="mm", res=900, pointsize=7)

ggplot(Percentages, aes(x=year, y=HSC)) + 
  geom_line(data = Percentages, aes(x=year, y=HSC_SMO_per, col = "Higher-Grade\n Service Class\n and Owners", linetype = "Higher-Grade\n Service Class\n and Owners"), size=0.5) +
  geom_line(data = Percentages, aes(x=year, y=LSC, col = "Lower-Grade\n Service Class", linetype = "Lower-Grade\n Service Class"), size=0.5) +
  geom_line(data = Percentages, aes(x=year, y=SW, col = "Skilled\n Workers", linetype = "Skilled\n Workers"), size=0.5) +
  geom_line(data = Percentages, aes(x=year, y=UW, col = "Unskilled\n Workers", linetype = "Unskilled\n Workers"), size=0.5) +
  ylab("Proportion") + xlab("") + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.5)
  ) + ggtitle("") + theme(plot.title = element_text(hjust = 0.5)) +  theme(text = element_text(size = 7, family="sans")) + 
  coord_cartesian(ylim = c(0,100)) +
  theme(legend.position="bottom") + theme(legend.key.width = unit(1,"cm")) +
  scale_color_manual(name = "", values = c("black", "gray25", "gray50", "gray75"), labels = c("Higher-Grade \nService Class and Owners", "Lower-Grade \nService Class", "Skilled \nWorkers", "Unskilled \nWorkers")) +
  scale_linetype_manual(name = "", values = c("solid", "dotted", "dotdash", "dashed"), labels = c("Higher-Grade \nService Class and Owners", "Lower-Grade \nService Class", "Skilled \nWorkers", "Unskilled \nWorkers")) 

dev.off()


##################################################
#################### FIGURE 1 ####################
##################################################

jpeg(file="figure1.jpg", width=114, height=80, units="mm", res=900, pointsize=7)

p1 = ggplot(Dataset_CIS, aes(x=tertiary, y=primary, color = as.factor(aprobada*-1))) + 
  geom_point(size=0.5) + geom_line(aes(y = 50)) +  ggtitle("Education") +
  geom_line(aes(x = 50)) + labs (x="Support among tertiary\n education (%)", y="Support among primary education (%)") +
  scale_color_grey() + theme_bw(base_line_size = 0.5) + theme(legend.position = "none") + theme(plot.title = element_text(hjust = 0.5), 
                                                                                                panel.border = element_rect(colour = "black", fill=NA, size=0.5),
                                                                                                panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(text = element_text(size = 7, family="sans"))

p2 = ggplot(Dataset_CIS, aes(x=HSC_SMO_pref, y=unskilled_workers, color=as.factor(aprobada*-1))) + 
  geom_point(size=0.5) + geom_line(aes(y = 50)) + ggtitle("Occupation") +
  geom_line(aes(x = 50)) + labs (x="Support among higher-grade\n service class and owners (%)", y="Support among unskilled workers (%)") +
  scale_color_grey() + theme_bw(base_line_size = 0.5) + theme(legend.position = "none") + theme(plot.title = element_text(hjust = 0.5), 
                                                                                 panel.border = element_rect(colour = "black", fill=NA, size=0.5),
                                                                                 panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(text = element_text(size = 7, family="sans"))

grid.arrange(p1, p2, nrow = 1)

dev.off()


##################################################
#################### FIGURE 2 ####################
##################################################

pref_gap_education_abs <- na.omit(abs(pref_gap_education))
mean(pref_gap_education_abs, na.rm=T)

pref_gap_five_abs <- abs(pref_gap_five)
mean(pref_gap_five_abs, na.rm=T)


jpeg(file="figure2.jpg", width=114, height=80, units="mm", res=900, pointsize=7)

plot1 = ggplot(Dataset_CIS, aes(y=pref_gap_education_abs, x=familia)) + 
  geom_bar(position="dodge", stat = "summary", fun.y = "mean", fill = "gray") +
  ylab("Preference gap") + xlab("") + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.5)
            ) + ggtitle("Education") + theme(plot.title = element_text(hjust = 0.5)) +  theme(text = element_text(size = 7, family="sans")) + 
  theme(axis.title.y = element_text(size=7, margin = margin(t = 0, r = 10, b = 0, l = 0))) +
  theme(axis.title.x = element_text(size=7, margin = margin(t = 0, r = 0, b = 0, l = 0)), axis.text.x = element_text(angle = 45, hjust = 1)) + 
  scale_x_discrete(labels=c("Economy", "Political processes", "Security and\n foreign affairs", "Cultural issues")) +
  scale_y_continuous(labels = scales::number_format(accuracy = 1)) + 
  coord_cartesian(ylim = c(0,20))

plot2 = ggplot(Dataset_CIS, aes(y=pref_gap_five_abs, x=familia)) + 
  geom_bar(position="dodge",  stat = "summary", fun.y = "mean", fill = "gray") + 
  ylab("Preference gap") + xlab("") + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.5)
                     ) + ggtitle("Occupation") + theme(plot.title = element_text(hjust = 0.5)) + theme(text = element_text(size = 7, family="sans")) + 
  theme(axis.title.y = element_text(size=7, margin = margin(t = 0, r = 10, b = 0, l = 0))) +
  theme(axis.title.x = element_text(size=7, margin = margin(t = 0, r = 0, b = 0, l = 0)), axis.text.x = element_text(angle = 45, hjust = 1)) + 
  scale_x_discrete(labels=c("Economy", "Political processes", "Security and\n foreign affairs", "Cultural issues")) +
  scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
  coord_cartesian(ylim = c(0,20))

grid.arrange(plot1, plot2, ncol=2)

dev.off()


##################################################
################### TABLE A1 #####################
##################################################

model_1 <- glm(aprobada ~ tertiary, data = Dataset_CIS, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ secondary, data = Dataset_CIS, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ primary, data = Dataset_CIS, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ tertiary + secondary + primary, data = Dataset_CIS, family = binomial(link = "logit"))


##################################################
################### TABLE A2 #####################
##################################################

model_1 <- glm(aprobada ~ HSC_SMO_pref, data = Dataset_CIS, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ low_service_class, data = Dataset_CIS, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ skilled_workers, data = Dataset_CIS, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ unskilled_workers, data = Dataset_CIS, family = binomial(link = "logit"))

model_5 <- glm(aprobada ~ HSC_SMO_pref + low_service_class + skilled_workers + unskilled_workers, data = Dataset_CIS, family = binomial(link = "logit"))


##################################################
################### FIGURE 3 #####################
##################################################

quantile(tertiary, c(.25, .75), na.rm=T) 
quantile(secondary, c(.25, .75), na.rm=T) 
quantile(primary, c(.25, .75), na.rm=T) 

fullmodel <- glm(aprobada ~ tertiary + secondary + primary, data = Dataset_CIS, family = binomial(link = "logit"))
dat01 <- ggpredict(fullmodel, terms="tertiary[0:100 by=.01]")
dat02 <- ggpredict(fullmodel, terms="secondary[0:100 by=.01]")
dat03 <- ggpredict(fullmodel, terms="primary[0:100 by=.01]")

quantile(HSC_SMO_pref, c(.25, .75), na.rm=T) 
quantile(low_service_class, c(.25, .75), na.rm=T) 
quantile(skilled_workers, c(.25, .75), na.rm=T) 
quantile(unskilled_workers, c(.25, .75), na.rm=T) 

fullmodel2 <- glm(aprobada ~ HSC_SMO_pref + low_service_class + skilled_workers + unskilled_workers, data = Dataset_CIS, family = binomial(link = "logit"))
dat1 <- ggpredict(fullmodel2, terms="HSC_SMO_pref[0:100 by=.01]")
dat2 <- ggpredict(fullmodel2, terms="low_service_class[0:100 by=.01]")
dat3 <- ggpredict(fullmodel2, terms="skilled_workers[0:100 by=.01]")
dat4 <- ggpredict(fullmodel2, terms="unskilled_workers[0:100 by=.01]")


jpeg(file="figure3.jpg", width=114, height=80, units="mm", res=900, pointsize=7)
par(mfrow=c(1,2), mar=c(4, 5, 2, 1) + 0.1)

plot(c(0,100),c(0,120),type="n",axes=F,
     ylim=range(0:120),
     xlim=range(0:100),
     ylab="Probability of policy change",
     xlab="Proportion support",
     main="Education"
)
box()
labs <- seq(0, 100, by = 25)
axis(1, at = labs, labels = paste0(labs, "%"), cex=0.8)
axis(2, at = labs, labels = paste0(labs, "%"))
lines(dat01$x, dat01$predicted*100, col="black", lwd=2)
lines(dat01$x[dat01$x>37.6225 & dat01$x<82.625], dat01$predicted[dat01$x>37.6225 & dat01$x<82.625]*100, col="black", lwd=3)
lines(dat02$x, dat02$predicted*100, col="gray50", lwd=2)
lines(dat02$x[dat02$x>38.61562 & dat02$x<78.9225], dat02$predicted[dat02$x>38.61562 & dat02$x<78.9225]*100, col="gray50", lwd=3)
lines(dat03$x, dat03$predicted*100, col="gray75", lwd=2)
lines(dat03$x[dat03$x>35.315 & dat03$x<81.4475], dat03$predicted[dat03$x>35.3150 & dat03$x<81.4475]*100, col="gray75", lwd=3)
legend(x="topright", cex=0.7, c("Tertiary","Secondary","Primary"), col=c("black", "gray50","gray75"), lwd=3)


plot(c(0,100),c(0,120),type="n",axes=F,
     ylim=range(0:120),
     xlim=range(0:100),
     ylab="Probability of policy change",
     xlab="Proportion support",
     main="Occupation"
)
box()
labs <- seq(0, 100, by = 25)
axis(1, at = labs, labels = paste0(labs, "%"), cex=0.8)
axis(2, at = labs, labels = paste0(labs, "%"))
lines(dat1$x, dat1$predicted*100, col="black", lwd=2)
lines(dat1$x[dat1$x>42.7950 & dat1$x<78.4325], dat1$predicted[dat1$x>42.7950 & dat1$x<78.4325]*100, col="black", lwd=3)
lines(dat2$x, dat2$predicted*100, col="gray25", lwd=2)
lines(dat2$x[dat2$x>42.3200 & dat2$x<81.2725], dat2$predicted[dat2$x>42.3200 & dat2$x<81.2725]*100, col="gray25", lwd=3)
lines(dat3$x, dat3$predicted*100, col="gray50", lwd=2)
lines(dat3$x[dat3$x>41.3550 & dat3$x<79.6275], dat3$predicted[dat3$x>41.3550 & dat3$x<79.6275]*100, col="gray50", lwd=3)
lines(dat4$x, dat4$predicted*100, col="gray75", lwd=2)
lines(dat4$x[dat4$x>41.9725 & dat4$x<80.9975], dat4$predicted[dat4$x>41.9725 & dat4$x<80.9975]*100, col="gray75", lwd=3)
legend(x="topright", cex=0.7, c("Higher-grade service & owners","Lower-grade service","Skilled workers","Unskilled workers"), lwd=3, col=c("black", "gray25", "gray50","gray75"))

dev.off()


##################################################
################### FIGURE A3 ####################
##################################################

fullmodel <- glm(aprobada ~ tertiary + secondary + primary, data = Dataset_CIS, family = binomial(link = "logit"))
dat01 <- ggpredict(fullmodel, terms="tertiary[0:100 by=.01]")
dat02 <- ggpredict(fullmodel, terms="secondary[0:100 by=.01]")
dat03 <- ggpredict(fullmodel, terms="primary[0:100 by=.01]")

fullmodel2 <- glm(aprobada ~ HSC_SMO_pref + low_service_class + skilled_workers + unskilled_workers, data = Dataset_CIS, family = binomial(link = "logit"))
dat1 <- ggpredict(fullmodel2, terms="HSC_SMO_pref[0:100 by=.01]")
dat2 <- ggpredict(fullmodel2, terms="low_service_class[0:100 by=.01]")
dat3 <- ggpredict(fullmodel2, terms="skilled_workers[0:100 by=.01]")
dat4 <- ggpredict(fullmodel2, terms="unskilled_workers[0:100 by=.01]")


jpeg(file="figurea3.jpg", width=6.5, height=4.5, units="in", res=900, pointsize=9)
par(mfrow=c(1,2), mar=c(4, 5, 2, 0) + 0.1)

plot(c(0,100),c(0,120),type="n",axes=F,
     ylim=range(0:120),
     xlim=range(0:100),
     ylab="Probability of policy change",
     xlab="Proportion support",
     main="Education", cex.main=1.5, font.main=1
)
box()
labs <- seq(0, 100, by = 25)
axis(1, at = labs, labels = paste0(labs, "%"), cex=0.8)
axis(2, at = labs, labels = paste0(labs, "%"))
lines(dat01$x, dat01$predicted*100, col="black", lwd=3)
lines(dat01$x, dat01$conf.low*100, col="black", lwd=1, lty="dotdash")
lines(dat01$x, dat01$conf.high*100, col="black", lwd=1, lty="dotdash")
lines(dat02$x, dat02$predicted*100, col="gray50", lwd=3)
lines(dat02$x, dat02$conf.low*100, col="gray50", lwd=1, lty="dotdash")
lines(dat02$x, dat02$conf.high*100, col="gray50", lwd=1, lty="dotdash")
lines(dat03$x, dat03$predicted*100, col="gray75", lwd=3)
lines(dat03$x, dat03$conf.low*100, col="gray75", lwd=1, lty="dotdash")
lines(dat03$x, dat03$conf.high*100, col="gray75", lwd=1, lty="dotdash")
legend(x="topright", cex=0.9, c("Tertiary","Secondary","Primary"), col=c("black", "gray50","gray75"), lwd=3)


plot(c(0,100),c(0,120),type="n",axes=F,
     ylim=range(0:120),
     xlim=range(0:100),
     ylab="Probability of policy change",
     xlab="Proportion support",
     main="Occupation", cex.main=1.5, font.main=1
)
box()
labs <- seq(0, 100, by = 25)
axis(1, at = labs, labels = paste0(labs, "%"), cex=0.8)
axis(2, at = labs, labels = paste0(labs, "%"))
lines(dat1$x, dat1$predicted*100, col="black", lwd=3)
lines(dat1$x, dat1$conf.low*100, col="black", lwd=1, lty="dotdash")
lines(dat1$x, dat1$conf.high*100, col="black", lwd=1, lty="dotdash")
lines(dat2$x, dat2$predicted*100, col="gray25", lwd=3)
lines(dat2$x, dat2$conf.low*100, col="gray25", lwd=1, lty="dotdash")
lines(dat2$x, dat2$conf.high*100, col="gray25", lwd=1, lty="dotdash")
lines(dat3$x, dat3$predicted*100, col="gray50", lwd=3)
lines(dat3$x, dat3$conf.low*100, col="gray50", lwd=1, lty="dotdash")
lines(dat3$x, dat3$conf.high*100, col="gray50", lwd=1, lty="dotdash")
lines(dat4$x, dat4$predicted*100, col="gray75", lwd=3)
lines(dat4$x, dat4$conf.low*100, col="gray75", lwd=1, lty="dotdash")
lines(dat4$x, dat4$conf.high*100, col="gray75", lwd=1, lty="dotdash")
legend(x="topright", cex=0.9, c("Higher-grade service & owners","Lower-grade service","Skilled workers","Unskilled workers"), lwd=3, col=c("black", "gray25", "gray50","gray75"))

dev.off()



##################################################
################### TABLE 1 ######################
##################################################

pop_tertiary = Dataset_CIS[which(tertiary>=75),]
mean(pop_tertiary$aprobada, na.rm=T)
nrow(pop_tertiary)

pop_primary = Dataset_CIS[which(primary>=75),]
mean(pop_primary$aprobada, na.rm=T)
nrow(pop_primary)

morepop_tertiary = Dataset_CIS[which(tertiary>=75 & pref_gap_education>5),]
mean(morepop_tertiary$aprobada, na.rm=T)
nrow(morepop_tertiary)

morepop_primary = Dataset_CIS[which(primary>=75 & pref_gap_education < -5),]
mean(morepop_primary$aprobada, na.rm=T)
nrow(morepop_primary)

pop_HSC_SMO = Dataset_CIS[which(HSC_SMO_pref>=75),]
mean(pop_HSC_SMO$aprobada, na.rm=T)
nrow(pop_HSC_SMO)

pop_unskilled_workers = Dataset_CIS[which(unskilled_workers>=75),]
mean(pop_unskilled_workers$aprobada, na.rm=T)
nrow(pop_unskilled_workers)

morepop_HSC_SMO = Dataset_CIS[which(HSC_SMO_pref>=75 & pref_gap_five>5),]
mean(morepop_HSC_SMO$aprobada, na.rm=T)
nrow(morepop_HSC_SMO)

morepop_unskilled_workers = Dataset_CIS[which(unskilled_workers>=75 & pref_gap_five < -5),]
mean(morepop_unskilled_workers$aprobada, na.rm=T)
nrow(morepop_unskilled_workers)


##################################################
################### TABLE A3 #####################
##################################################

disagree_educ5 <- subset(Dataset_CIS, pref_gap_education_abs>=5)
disagree_educ10 <- subset(Dataset_CIS, pref_gap_education_abs>=10)
disagree_ocup5 <- subset(Dataset_CIS, pref_gap_five_abs>=5)
disagree_ocup10 <- subset(Dataset_CIS, pref_gap_five_abs>=10)


model_1 <- glm(aprobada ~ pref_gap_education, data = Dataset_CIS, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ pref_gap_education, data = disagree_educ5, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ pref_gap_education, data = disagree_educ10, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ pref_gap_five, data = Dataset_CIS, family = binomial(link = "logit"))

model_5 <- glm(aprobada ~ pref_gap_five, data = disagree_ocup5, family = binomial(link = "logit"))

model_6 <- glm(aprobada ~ pref_gap_five, data = disagree_ocup10, family = binomial(link = "logit"))



##################################################
#################### FIGURE 4 ####################
##################################################

jpeg(file="figure4.jpg", width=114, height=80, units="mm", res=900, pointsize=7)

plot1 = effect_plot(model_1, pred = pref_gap_education, interval = TRUE, 
                    y.label = "Predicted probablity of policy change", x.label="Preference gap", plot.points = FALSE, main.title="Education") +
                    scale_y_continuous(labels = scales::percent, limits=c(0,1)) + 
                    drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
                    theme_bw(base_line_size = 0.5)  + theme(text = element_text(size = 7, family="sans")) + 
                    theme(plot.title = element_text(hjust = 0.5), panel.border = element_rect(colour = "black", fill=NA, size=0.25), 
                          panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
                    theme(axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0)))
  
plot2 = effect_plot(model_4, pred = pref_gap_five, interval = TRUE, 
                    y.label = "Predicted probablity of policy change", x.label="Preference gap", plot.points = FALSE, main.title="Occupation") + 
                    scale_y_continuous(labels = scales::percent, limits=c(0,1)) + 
                    drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
                    theme_bw(base_line_size = 0.5)  + theme(text = element_text(size = 7, family="sans")) + 
                    theme(plot.title = element_text(hjust = 0.5), panel.border = element_rect(colour = "black", fill=NA, size=0.25), 
                          panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
                    theme(axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0)))
  
grid.arrange(plot1, plot2, ncol=2)

dev.off()



##################################################
################### TABLE A4 #####################
##################################################

economy <- subset(Dataset_CIS, familia=="Economic Issues")
politics <- subset(Dataset_CIS, familia=="Politics")
security <- subset(Dataset_CIS, familia=="Security and Foreign Affairs")
social <- subset(Dataset_CIS, familia=="Cultural Issues")

model_1 <- glm(aprobada ~ pref_gap_education, data = economy, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ pref_gap_education, data = politics, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ pref_gap_education, data = security, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ pref_gap_education, data = social, family = binomial(link = "logit"))


##################################################
#################### FIGURE 5 ####################
##################################################

jpeg(file="figure5.jpg", width=114, height=176, units="mm", res=900, pointsize=7)

plot1 = effect_plot(model_1, pred = pref_gap_education, interval = TRUE, 
                    y.label = "Predicted probablity of policy change", x.label="Preference gap", plot.points = FALSE, main.title="Economy") +
  scale_y_continuous(labels = scales::percent, limits=c(0,1)) + 
  drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
  theme_bw(base_line_size = 0.5)  + theme(text = element_text(size = 7, family="sans")) +
  theme(plot.title = element_text(hjust = 0.5), panel.border = element_rect(colour = "black", fill=NA, size=0.25), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.title.y = element_text(size=7, margin = margin(t = 0, r = 5, b = 0, l = 0)))

plot2 = effect_plot(model_2, pred = pref_gap_education, interval = TRUE, 
                    y.label = "Predicted probablity of policy change", x.label="Preference gap", plot.points = FALSE, main.title="Political Processes") +
  scale_y_continuous(labels = scales::percent, limits=c(0,1)) + 
  drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
  theme_bw(base_line_size = 0.5)  + theme(text = element_text(size = 7, family="sans")) +
  theme(plot.title = element_text(hjust = 0.5), panel.border = element_rect(colour = "black", fill=NA, size=0.25), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.title.y = element_text(size = 7, margin = margin(t = 0, r = 5, b = 0, l = 0)))

plot3 = effect_plot(model_3, pred = pref_gap_education, interval = TRUE, 
                    y.label = "Predicted probablity of policy change", x.label="Preference gap", plot.points = FALSE, main.title="Security and Foreign Affairs") +
  scale_y_continuous(labels = scales::percent, limits=c(0,1)) + 
  drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
  theme_bw(base_line_size = 0.5)  + theme(text = element_text(size = 7, family="sans")) +
  theme(plot.title = element_text(hjust = 0.5), panel.border = element_rect(colour = "black", fill=NA, size=0.25), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.title.y = element_text(size = 7, margin = margin(t = 0, r = 5, b = 0, l = 0)))

plot4 = effect_plot(model_4, pred = pref_gap_education, interval = TRUE, 
                    y.label = "Predicted probablity of policy change", x.label="Preference gap", plot.points = FALSE, main.title="Cultural Issues") +
  scale_y_continuous(labels = scales::percent, limits=c(0,1)) + 
  drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
  theme_bw(base_line_size = 0.5)  + theme(text = element_text(size = 7, family="sans")) +
  theme(plot.title = element_text(hjust = 0.5), panel.border = element_rect(colour = "black", fill=NA, size=0.25), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.title.y = element_text(size = 7, margin = margin(t = 0, r = 5, b = 0, l = 0)))

grid.arrange(plot1, plot2, plot3, plot4, ncol=2)

dev.off()



##################################################
################### TABLE A5 #####################
##################################################

model_1 <- glm(aprobada ~ pref_gap_education + gov_party + pref_gap_education*gov_party, data = Dataset_CIS, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ pref_gap_five + gov_party + pref_gap_five*gov_party, data = Dataset_CIS, family = binomial(link = "logit"))


##################################################
#################### FIGURE 6 ####################
##################################################

jpeg(file="figure6.jpg", width=114, height=80, units="mm", res=900, pointsize=7)

plot1 = interplot(m = model_1, var1 = "pref_gap_education", var2 = "gov_party") + xlab("Government ideology") +
  ylab("Predicted effect of preference gap") +
  scale_x_continuous(breaks=c(0, 1), labels=c("Left", "Right"), limits = c(-.5, 1.5)) +
  scale_y_continuous(limits = c(-0.2, 0.2)) +
  drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
  geom_hline(yintercept=0, linetype="dotted") + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.25)) + 
  ggtitle("Education") + theme(plot.title = element_text(hjust = 0.5)) +  
  theme(axis.title.y = element_text(size = 10, margin = margin(t = 0, r = 5, b = 0, l = 0)))


plot2 = interplot(m = model_2, var1 = "pref_gap_five", var2 = "gov_party") + xlab("Government ideology") +
  ylab("Predicted effect of preference gap") +
  scale_x_continuous(breaks=c(0, 1), labels=c("Left", "Right"), limits = c(-.5, 1.5)) +
  scale_y_continuous(limits = c(-0.2, 0.2)) +
  drop_gridlines(x = TRUE, y = TRUE, minor.only = FALSE) +
  geom_hline(yintercept=0, linetype="dotted") + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.25)) + 
  ggtitle("Occupation") + theme(plot.title = element_text(hjust = 0.5)) +  
  theme(axis.title.y = element_text(size = 10, margin = margin(t = 0, r = 5, b = 0, l = 0)))

grid.arrange(plot1, plot2, ncol=2)

dev.off()



##################################################
################### TABLE A6 #####################
##################################################

allincome <- subset(Dataset_CIS, per_ing_5>=0)

model_1 <- glm(aprobada ~ per_ing_5, data = allincome, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ per_ing_3, data = allincome, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ per_ing_1, data = allincome, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ per_ing_5 + per_ing_3 + per_ing_1, data = allincome, family = binomial(link = "logit"))

model_5 <- glm(aprobada ~ pref_gap_income2, data = allincome, family = binomial(link = "logit"))


##################################################
################### TABLE A7 #####################
##################################################

noindirect <- subset(Dataset_CIS, indirect_question==0)
since1993 <- subset(Dataset_CIS, year>1993)

model_1 <- glm(aprobada ~ pref_gap_education, data = noindirect, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ pref_gap_five, data = noindirect, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ pref_gap_education, data = since1993, family = binomial(link = "logit"))


##################################################
################### TABLE A8 #####################
##################################################

model_1 <- glm(aprobada ~ pref_gap_education + realgdpgr + unemp + inflation + rae_leg + gov_party, data = Dataset_CIS, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ pref_gap_education + realgdpgr + unemp + inflation + rae_leg + gov_party, data = disagree_educ5, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ pref_gap_education + realgdpgr + unemp + inflation + rae_leg + gov_party, data = disagree_educ10, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ pref_gap_five + realgdpgr + unemp + inflation + rae_leg + gov_party, data = Dataset_CIS, family = binomial(link = "logit"))

model_5 <- glm(aprobada ~ pref_gap_five + realgdpgr + unemp + inflation + rae_leg + gov_party, data = disagree_ocup5, family = binomial(link = "logit"))

model_6 <- glm(aprobada ~ pref_gap_five + realgdpgr + unemp + inflation + rae_leg + gov_party, data = disagree_ocup10, family = binomial(link = "logit"))



##################################################
################### TABLE A9 #####################
##################################################

model_1 <- glm(aprobada ~ pref_gap_five, data = economy, family = binomial(link = "logit"))

model_2 <- glm(aprobada ~ pref_gap_five, data = politics, family = binomial(link = "logit"))

model_3 <- glm(aprobada ~ pref_gap_five, data = security, family = binomial(link = "logit"))

model_4 <- glm(aprobada ~ pref_gap_five, data = social, family = binomial(link = "logit"))

